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STATEMENT  OF  PROBLEM; 


Let  x:  a  =  x^<X2< - =  b 


=  d 


y  ;  c  =  y^<  y2<l - <yjj 

and  n  =  [a,bj  x  [c,d]  . 

Given  (probably  noisy)  data 

flj  =  fUi.yj)  .  tij,  1il,j£N 

let  0<^<1,  and  let  vf--'yO 
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minimize 
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dxdy 


over  appropriate  smooth  functions  g. 


SUMMARY  OP  RESULTS; 


The  investigation  of  the  problem  has  included  the  work  of 

[1|2  ,  • 

As  in  the  one  dimensional  case,  Chui  points  out  that  it  is 
sufficient  to  minimize  the  auxiliary  expression 


b  d  m  m 

f  f  7  k  dxdy 

L  L  k=0  ax^ay®”^ 


(1) 


over  appropriate  smooth  functions  g. 

As  analyzed  inC6»7l  a  proper  abstract  setting  for  a  natural 

2-dimensional  setting  of  this  optimal  problem  and  results  is 

provided  by  the  Beppo  Levi  space  X  =  BL™(R^)  of  order  m  over 
2 

R  (m  an  integer  greater  than  or  equal  1).  Defined  as 
X  =  {g€D';  a‘^g€L2  for  \oc\  =  m},. 

where  <?<-*  - and|o<|  =  o(^  +o|  + - +64, 

X  is  thus  simply  the  vector  space  of  all  the  (Schwartz) 


(2) 


distributions  (i.e.  continuous  linear  functionals  on  the  vector 

space  of  infinitely  differentiable  functions  with  compact 

support  in  R  ,  provided  with  the  canonical  Schwartz  topology) 

for  which  all  the  partial  derivatives  (in  the  distributional 

2 

sense)  of  total  order  m  are  square  integrable  in  R  •  X  is 
naturally  equipped  with  the  semi-inner  product  C*t*)m  correspond 
ing  to  the  rotation  invariant  seminorm 


the  Kernel  of  I  *  known  to  be  simply  the  vector  space 

P  =  ^  of  dimension 

m— I 


2 

of  all  polynomials  over  R  of  (total)  degree  m  -  1,  It  should 

2 

be  noted  that  igj^j^  may  be  physically  interpreted  (at  least  if 
m  =  2  and  under  some  simplifying  assumptions)  as  the  bending 
energy  of  a  thin  plate  of  infinite  extent,  g  denoting  the 
deflection  normal  to  the  rest  position.  As  they  define  the 
equilibrium  position  of  infinite  extent  that  deforms  in  bending 
only  (under  deflections  specified  at  a  number  of  independent 
points),  the  solution  of  the  optimed  interpolation  problem  can 
be  appropriately  termed  surface  splines. 

Originally  introduced  for  interpolating  wing  deflections 
and  computing  slopes  for  aeroelastic  calculations  by  (51,  this 
ingenious  device  proves  most  interesting  to  analyze  mathemati¬ 
cally;  in  this  connection,  various  deep  results  have  been  ob¬ 
tained  inC2l.  However  inc63  a  more  constructive  approach  is 
taken,  where  a  prominent  role  is  played  by  representation 
formulas  in  function  and  distribution  spaces,  these  comple¬ 
mentary  results  being  obtained  by  resorting  to  such  basic 


mathematical  tools  as  convolutions  and  Fourier  transforms  of 
distributions.  For  a  more  concrete  presentation  stressing  the 
significant  properties  of  the  optimal  interpolation  process  of 
surface  spline  interpolation  at  the  algorithmic  level  we  present 
the  approach  of  Meinguet  [ 7  J  . 

Meinguet  formulates  the  optimal  problem  as  follows: 

Let  there  be  given: 

A  finite  set  A  =  of  distinct  points  of  R  contain¬ 

ing  a  p  -unisolvent  subset,  by  which  we  mean  a  set  B  =  Caj)j^j 
of  M  points  of  A,  M  being  defined  by  (4),  such  that  there  exists 
a  unique  pcP satisfying  the  interpolating  conditions 

P(aj)  ,  V  3  e  J  (5) 

for  any  prescribed  real  scalars  ©4^ ,  V  j  eJ. 

il 

A  set  of  real  scalars  or  equivalently  provided 

that  m  >  S  ,  the  linear  variety : 


V  =  £gcX:  g(a^)^  ,  V  i€l)  ;  (6) 

whenever  s  f(aj^),  Vi€-I,  where  f  denotes  a  function 
defined  (at  least)  on  A,  V  can  be  interpreted  as  the  set  of 
X-interpolants  of  f  on  A.  Thus  our  problem  is  to  find  h£V 
such  that 


(7) 


By  virtue  of  the  p  -unisolvence  of  the  subset  B  of  the  given 

2 

set  A  of  interpolation  points  in  R  ,  there  exists  in  P=Pni_i 
a  unique  basis  (Pj)i  j  that  is  dual  to  the  set  of  shifted 
dirac  measures  (in  the  sense  that  pj^(aj)  =  i»3<^J 

where  is  the  Kronecker  symbol).  For  every  g-eX,  ra>1 ,  the 


P-interpolant  P  of  g  on  B  is  given  by 

s 


‘g 


(8) 


owing  to  this  definition,  the  mapping  P:  X-»X  is  a  linear 
projection  of  X  with  range  P=P„,.^ 


and  Kernel 


(9) 


=  [ge  X:  g(a.)  =  0,  V  j  €  J}, 


so  that 


*  =  Vl  ® 

equipped  with  the  seminorm  |  ^  Hilbert  space. 

In  view  of  the  direct  sum  decomposition  (10),  the  restrict¬ 
ion  of  the  associated  projector  I- P  to  the  linear  variety  VCX 
defined  by  (6)  is  clearly  an  injection.  Therefore,  finding  heV 
such  that  (7)  holds  amounts  strictly  to  finding  an  element  w  of 
minimal  norm  |  •  l  in  the  image  of  V  under  I  — P,  which  is  the 

linear  variety, 

W  =  (g6X^:  g(aj^)  =o^,VkaK},  (11) 

where 

Kil- J  =  lk6  N:  1  ^k^N=  card  (K)}  (12) 

for  definiteness,  and 

The  unique  solution  h  of  the  problem  is  given  by 
h  =  w  + 


3^3 


(U) 


depend  conjiinuously  on  the  data  (<^i)ifex  *  where 

«  -  &  ■  ('5) 

denotes  the  Prechet:-Riesz  representer  of  the  a  shifted 
Di^ac  measure  c(_  ,  YL  real  coefficients  satisfying  the  Cramer 

system  of  linear  equations 

^  K(a.,a^)\  =o(^  ,  1^k^N.  (16) 

fel 


in  fact  K_  involves  no  functions  more  complicated  than  logarithnxs 

^  2 
and  is  easily  coded.  The  set  {K^:Vx€R  }  is  the  so  called 

reproducing  kernel  of  the  Hilbert  function  space  X^;  it  can  be 

regarded  equivalently  as  the  real-valued  continuous  function 

(x,y)-»K(x,y)  =  K^(y)  on  R^xR^. 

Two  dimensional  interpolation  by  radial  basis  functions. 

The  ’’thin  plate  spline”  is  among  the  radial  basis  functions  for 


surface  interpolation  which  were  found  to  be  most  successful  ^4] . 
The  drawback  in  using  these  basis  functions  for  interpolation  is 
the  need  to  solve  large  full  systems  of  linear  equations  which 
become  very  ill-conditioned  as  their  order,  namely  the  number  of 
data  points,  increases.  However,  the  abov^e  basis  function  share 
the  property  that 


^  wCr)- 


as  r  increases 


(r^  =  +  y^) 


and  41  w(r)->co  as  r-»0 


where  a  is  the  m-iterated  Laplacian,A= 


(17) 


This  property  of  the  basis  functions  enables  us  to  apply  a  pro¬ 
cedure  by  which  most  of  the  interpolation  equations  become 
diagonally  dominant. 

As  pointed  out  in  [  33 ,  the  major  application  of  the  above 
radial  basis  functions  is  for  interpolation  of  scattered  data. 
This  presentation  considers  data  given  on  a  square  grid,  where 
the  discrete  analogues  of  the  operators  a"*  are  well  known.  A 
method  for  defining  the  discrete  analogues  of  a”  for  general 
domain  is  being  investigated. 

Let  the  data  f^^^  =  f(x^,y^)  be  given  on  the  square  grid: 

(Xi,y.)  =  ((i-1)h,(j-1)h),  1^i,j^N  (18) 


The  system  of  interpolation  equations  becomes 
N  ,  . 


(19) 


where 


(Xi-Xk)2  +  (y.-yi)' 


We  would  like,  by  row  operations,  to  form  a  finite  difference 
approximation  to  the  iterated  Laplacian  ofw(r),  and  thus  to 
generate  diagonal  dominance  in  the  system  (19).  The  central 
difference  operator  ^  is  replaced  here  by  the  5-point 
difference  approximation  to  the  Laplacian 

^h^kl  *  fb-j-l  1  fv  1  J.1  fv  1  _1  “  4fi,T  • 


*  ^k+1,1  ^k,l+1 


^k-1,1  ^k,l-l 


Operating  with  (with  respect  to  the  indices  k  and  1)  on 


(19)  we  obtain 
i»  0  =  1 


By  properties  of  (17)  we  expect  thatA®w(r^j)  decreases  as  r^^ 

increases  and  that  the  diagonal  element  becomes  large.  Further¬ 
more  we  expect  to  achieve  diagonal  dominance  in  the  transformed 
rows ,  i . e .  , 


N 


ASw(r“)l 


m=1  ^  k,l  ^  N-m. 


The  remaining  equations  cannot,  in  general,  be  made  diagonally 
dominant.  Therefore,  even  if  the  major  part  of  the  system  (19) 
can  be  transformed  into  a  diagonally  dominant  system  to  be  solved 
by  iterations,  there  is  a  substantial  part  of  the  system  which 
must  be  solved  directly  by  each  iteration.  A  way  of  overcoming 
this  difficulty,  by  using  special  differencing  of  (19)  near  the 
"boundaries”  of  the  system,  is  presented. 


Iterative  construction  of  "thin  plate  splines','  TPS,  on 
a  rectangular  grid. 

We  consider  TPS  with  ra=2;  this  surface  corresponds  to  an 

infinte  extent  thin  plate  of  minimum  bending  energy  clamped  at 

the  data  points  [3].  For  this  case  the  basis  functions  are 

2 

fundamental  solutions  of  the  biharmonic  equation  Aw  =  0,  i.e., 

2 

w(r)  =  r  log  r 

augmented  by  the  monomials  1,x  and  y.  The  interpolation  equations 
for  data  given  on  a  square  grid  (18)  cure 
N 

X  a  ♦  *  =^1  -  % 


1^k,l:^N,  (20) 


with  the  constraints 


N 


N 


N 


'  Z  '  Z  '  °' 


(21) 


k,l=1 
.2. 


k,l=1  k,l=1 

Since  A^w(r..)  =cr(r..),  it  is  expected  that  the  difference 

-L  J  J 

0 

operator  is  appropriate  for  generating  diagonal  dominance. 

So  far  the  process  of  differencing  the  equations  is  limited 
to  the  interior  equations  with  3^k,l^N-2.  The  application  is 


extended  to  the  "boundary”  equations,  i.e.,  with  k  or  1  equal  to 

2 

1,2,N-1,N  by  defining  the  difference  operator  on  the 
boundary  points  of  the  domain  in  such  a  way  that  =  0  on  all 


grid  points  if  and  only  if  f  is  a  linear  grid  function  of  the 
form  f^^  =  a  +  bx^  +  cy^,  1^k,l^N. 


This  property  guarantees  that  the  polynomial  part  in  (20) 

2 

vanishes  with  the  application  of  to  (20),  and  cam  be  restored 


from  the  solution  oi  of  the  derived  system; 
N 


A2 


Kp-I 

L  i,j  =  1 


kl 


=  0,  16k, 16N. 


(22) 


To  obtain  the  right  form  of  we  use  the  discrete  analogue  of 
the  iterated  Green’s  formula: 


a(f,g)s^(f^gxx+2fxy8xy*V«yy^' 

The  null  space  of  the  functional  a(f,f)  is  the  space  of  linear 
polynomials.  We  define  a  discrete  analogue  a^^  of  a,  so  that  the 
null  space  of  a^(f,f)  is  the  space  of  linear  functions: 

l,k=1 


k=' 


L=1 

N-2  N 


where  (<J^f  )tv  =  V  ^ i  v 


^‘^y^^lk^^yS^lk  ’ 


(23) 


Hence,  the  desired  property  of  is  guaranteed  by  formally 

defining  via  the  summation  by  parts  of  (23),  namely  by  the 
identity  jj. 

ah(f,g)  =^gkl^^h^^kl  * 

klT=1 

The  following  Lemmas  are  proved  int3l  : 

Lemma  1.  Let  the  difference  operator  be  defined  by  (23)  and 

2 

Then  Aj^f  =  0  if  and  only  if  f  is  a  grid  function  of  the 
form  f^j^  =  a  +  'bxj^+  cy^  ,  1  ^k,l^N. 


Lemma  2.  The  matrix  A  is  symmetric,  nonnegative  definite  of 
2 

rank  N  -3,  where  A  is  a  matrix  representation  of  the  difference 


operator  on  grid  functions,  regarded  as  vectors  in  R 

As  direct  conclusions  from  Lemmas  1  and  2,  we  have: 


N 

Corollary  1.  21  ’^kl*^^  ^^k  ^^1^  ~  ^  a,b,c  if 

k ,  1 "  1 
2 

and  only  if  w  for  some  grid  function  f. 


Corollary  2.  The oC-component  of  the  solution  of  (20)  and  (21) 
is  a  grid  function  of  the  form 

for  some  grid  function 

The  system  (22)  obtained  by  the  application  of to  (20)  can 
also  be  written  as 

'6k.l6Nor  (25) 

X  y  J  —  I 

7\Ao<  =Af, 

2  2  /  kX 

where  A  is  the  N  xN  matrix  with  entries  w(r.t),  ordered  in 

j 

accordance  with  the  vector  form  of  the  grid  functions.  The  sys¬ 
tem  (25)is  singular,  by  Lemma  2.  More  precisely  we  have: 
Corollary  3.  The  matrix  ^A  of  the  system  (25)  is  of  rank  N  -3. 

Due  to  the  symmetry  ofx^,  the  system  obtained  from  (20)  by 

2 

the  substitution 

Aj  *  *”'k  *  ■  ^ki- 


1^;k,l^rN,  (26) 


corresponds  to  the  representation  of  the  solution  in  terms  of  the 
new  basis  functions 

B^jCx.y)  '[A^wCr  _)]  14IJ6H, 

which  are  bell-shaped  for  3  6i,j^N-2.  This  system,  of  functions 

2 

is  of  dimension  N  -3,  by  corollary  1. 

In  particular: 

N  N  N 

Z!  Bi.Cx.y)  =  Z  x.B. .(x,y)  =  ^  y.B^.(x,y)  =  0. 
i7T=1^^  iTJ=^  ^  'iTT=l  ^  ^ 

The  solution  of  (25)  constitutes  a  three-dimensional  subspace 
2 

of  N  -  vectors  by  corollary  3.  In  order  to  obtain  that  solution  of 
(25)  which  satisfies  (21),  we  present  an  iterative  scheme  for  the 

solution  of  (25),  in  which  each  iterant  is  of  the  formAj^w,  and 


therefore  satisfies  (21)  in  view  of  corollary  1. 
The  iterative  scheme  is: 

=0  1^k,16N 

for  n  =  0, 1 ,2, - , 


^kl 


■  z  o<<“) 


w(rjb. 


1  £k,16  N 


(27) 

(28) 


]kl  •  (29) 

The  first  iterants  in  procedure  (29)  provide  smoothing 
solutions  for  noisy  data  {  f^^  . 

The  development  throughout  this  presentation  has  been 
restricted  to  squaire  grids.  However,  the  application  of  these 
ideas  to  triangular  grids  of  general  form  is  now  being  investi¬ 
gated.  In  addition  we  plan  to  continue  the  investigation  of  the 
smoothing  technique. 
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